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projected Gross-Pitaevskii equation to calculate correlation functions in 3D simulations of 
q-( varying trap elongation and temperature. 
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1. Introduction 

Macroscopic phase coherence is one of the key properties of a Bose-Einstein condensate 
(BEC), and has been a topic of longstanding interest. Almost as soon as BECs were created 
in the lab, this astounding property was demonstrated and investigated with respect to both 
spatial and temporal phase coherence El l2l El IH. Understanding phase evolution is key to 
the development of applications such as interferometry using condensates, and a great deal of 
work has been done in this area, see for example [El [6l • In addition, probing dynamic 
correlation functions by Bragg scattering is by now a mature technology jH HI. While 
much progress has been made, fundamental questions still remain over the nature of phase 
fluctuations with regard to the role of geometry and finite-temperature interactions ifTOlfTTI . 

The inclusion of a harmonic trapping potential adds complexity to the theoretical 
description of the gas, and the geometry of the trap dramatically alters the phase coherence 
properties ifTOl ITTL fT2"l| . In three dimensions, phase coherence can be spatially extensive 
and long-lived, even at finite temperatures iTOl ; in a tightly confined trap within the one 
dimensional limit, however, long range coherence is lost over the length of the condensate, 
and quasicondensates of length shorter than the system form. In this regime, theory has been 
developed lfT4l[T5l to calculate the coherence time of the system, and has seen good agreement 
with experiment f\A\. 

The cross-over regime between the ID and 3D limits is significantly more complex at 
finite temperature, and formalism is required that can accurately describe the full-dimensional 
trapped cloud, and capture the interactions in the system exactly. One appropriate class 
of formalisms are the so-called c-field methods, which have already had much success in 
calculating correlations within the finite temperature BEC |[T2"l[T3l[T6ll . 

In this work, we study the transition region between the ID and 3D limits. We make 
analytical calculations in the two extremes, and then develop formalism with the projected 
Gross Pitaevskii equation (PGPE) [171 [T8l to numerically study the cross-over region and 
compare with our analytic results. 

This paper is organised as follows. We begin in Section [2] with our system and correlation 
definitions. Bogoliubov analysis for the quasi ID and spherical limits is undertaken in Section 
[3l In Section |4] we outline the PGPE formalism we will use to make numerical studies of the 
finite temperature Bose gas. In Section [5] we give the results of our simulations for a variety 
of system symmetries and temperatures. We conclude in Section [6] 

2. System and Correlation Functions 

Our system is described by the second quantized Hamiltonian 




(1) 



where 



2m 



V 2 + H(r) 



(2) 
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is the single particle Hamiltonian, ^(r) is the quantum Bose field operator, and U = 
Aixh 2 a/m is the interaction strength, with a the s-wave scattering length. The trap potential is 
given as 

7 (r) = ~m (u 2 x x 2 + u 2 y y 2 + u 2 z z 2 ) , (3) 

where {cu x , u y ,co z } are the trap frequencies. 

In this work we will study the phase correlations within a system with respect to time. 
The phase correlation function between two points r 1; r 2 and two times ti, t 2 is given as 

g(n,t 1 ,r 2 ,t 2 ) = (^(n.tO^fo.fe)), (4) 

where the average is taken over realizations of field configurations. We also define a spatially 
averaged autocorrelation function, 

9n(t u t 2 ) = ^2g(r,t 1 ,T,t 2 ), (5) 

r 

where the points r are chosen to lie in the interior of the cloud where the density is high. 
This autocorrelation function gives the phase correlation of a single point in the system with 
respect to time, and we indicate this with the subscript '11'. We also define the two point 
autocorrelation function 

9i2(t u t 2 )= ^ 9( r uti, r 2,h), (6) 

<ri,r2> 

where both r x and r 2 are taken at interior points of system 1 and system 2, where the respective 
densities are high. This intercloud autocorrelation shows the phase correlation between two 
different systems with respect to time, and we signify this with the subscript '12'. 
We also define the normalised versions of Eqs. (I4I5I6I) as 

G[ti, ti, r 2 , t 2 ) = == (7) 

y/g{n,ti, ri, h)g(r 2l h, r 2 , t 2 ) 

Gix = ^= l ^^ ri,tl ' ri,t2 " > (8) 

n 

Lri2 — (?) 

n 

3. Bogoliubov Analysis 

In this section, we solve the Bogoliubov equations to find the decay of the autocorrelation 
within a cloud with respect to time. This is done in the quasi- ID and spherical limits. These 
analytical results will form the basis from which we will compare our numerical results. 

3.1. Quasi- ID Systems 

For a very elongated gas, the trap potential may be written 

V el {x, y, z) = l -m [co 2 ± (x 2 + y 2 ) + u 2 z z 2 ] . (10) 
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For later use, we define the asymmetry parameter 



A - — 



(11) 



and the oscillator lengths a± = (h/mu±) 1 ^ 2 and a z = a±/X 1 ^ 2 . We assume first A 1, so 
that the excitation of normal modes in the perpendicular direction can be neglected. For the 
case we consider, it is assumed that the interaction energy is high enough that the condensate 
occupies several trap modes, which enables us in the analytical calculations to treat the 
condensate wave function in the Thomas-Fermi approximation [fT9ll . 



n(x,y,z) 



9 



1 - 



x 2 + y 2 
R 2 



z 



(12) 



where 



R 2 
Z 2 



2fl 



2(JL 

moo 2 



2 ■ 



and finally 



ftw z fl5N a y /5 



V 



A 2 GU 



(13) 



(14) 



The calculation of autocorrelation functions for a Bose gas in the Bogoliubov 
approximation was outlined in Refs. [[8l [20l [2TJ . The correlation function is written 



g n (r,t,r',0) = y/n(r)n(r') exp 



-2 F M,r',0) 



where 



F(r,t,r',O) = (|0(r,t)-0(r',O)| 2 ) 
and the phase operator is defined as 



n r 



e^cij + h.c. 



(15) 



(16) 



(17) 



where finally /. (r) is a solution to the Bogoliubov equations in number-phase form ll22ll . 

[L-2n(r)]/r(r)=u^+(r), 

[L + gn(r)]f+{T) = ^/r(r), (18) 
and the linear operator L is defined as 

h 2 



-^-V 2 + V{r)-fi + 2gn(r). 
2m 



(19) 



The number-phase representation functions ff used here are related to the more familiar 
Bogoliubov amplitudes Uj and Vj through ft{ r ) — l u j( r ) i t 'j( r )]/2, where we have used the 
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sign and phase conventions of Ref. [fT9ll . The solution for the longitudinal modes, using the 
Thomas-Fermi approximation for the condensate density, is Il20ll23l 

-2gn(r) 1 (j + 2)(2j + 3)1 1/2 



Uj 8irR 2 Z j + 1 



uj =u z \J , (20) 

where j = 1,2,.. ., and P'- a ' h \x) is a Jacobi polynomial [|24l|. Computing the temporal 
correlation at the origin, r = r' = 0, for a temperature T assumed to be much larger than the 
longitudinal trap frequency u z , we obtain 

tiR l 7 + 1 a;,- a;,- 

j even 

It is found that at short times, t <C u^ 1 , F is proportional to the square of t, i.e. gn oc e~* 2 . 
For moderately short times, t ~ w^ 1 , F goes linear in t and the result for the correlation 
function is 

g(0,t,0,t)(xe- r ^ t , (22) 

where 

r _ 16k k B T n5Na\ 2/5 
15 tuu z \ \ 2 a z J 

For times t > u;" 1 , the correlation function will exhibit a quasiperiodic oscillatory behaviour. 
3.2. 3D Systems 

To contrast, we now treat analytically the case of a spherical condensate, with u± = u z = u. 
The excitation spectrum is given by E3l 

An + 3 2g 



1/2 



w n = un(2n + 3), (24) 

where n are positive integers. Following the method in the previous section, the exponent of 
the correlation function then becomes 

Ak B T 2g ^ An + 3 



hu ttR 3 u ^ \n(2n + 3)1 

n=l L ' 

(2n + l) 



2 2 ™(n! 



1 — cosu; n t). (25) 



At very short times, t <C u z 1 , F is once again proportional to the square of t. However, at 
short times t ~ w^ 1 F is approximately equal to 

\ 8qksT jut 

F <.°*°'*) = im?VT' (26) 
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so in the spherical case, the correlation goes as the square root of t for short times, and 
this distinction will be used in later sections to make comparisons between different trap 
symmetries. As in the previous case, when t > tu^ 1 , quasiperiodic behavior ensues. 

4. PGPE Formalism 

Motivated by the need to develop tractable theory to describe experimentally realistic systems, 
we now move to a system intermediary to the quasi ID and spherical 3D systems, that is 
to say, the elongated Bose gas. We develop the PGPE method to calculate temporal phase 
correlations in a three dimensional Bose gas at finite temperature. In this particular work, we 
note that it has a key advantage over Bogoliubov methods: the ease with which this method 
can accurately calculate temporal correlations for a wide variety of physical applications. We 
demonstrate this later in Section [5] by calculating results for differing trap aspect ratios and 
temperatures. 

In what follows, we give an overview of the PGPE method, and our sampling method 
for studying data. More details on the PGPE and c-field methods can be found in Refs. 



4.1. PGPE Formalism and Correlation Functions 

The key aspect of this formalism is that the Bose field operator is separated into two parts, 
dependent on the mean occupation of modes in the system. This splitting is written as 



where ipc is the coherent region c-field and t/>j is the incoherent field operator [17J. These 
fields are defined as the low and high energy projections of the full quantum field operator, 
separated by the energy e cut . In our theory this cutoff is implemented in terms of the 
harmonic oscillator eigenstates {ip n (r)} of the time-independent single particle Hamiltonian 
i.e. e n cp n (r) = H Q (p n (r), with e n the respective eigenvalues. 

In general, the applicability of the PGPE approach to describing the finite temperature 
gas relies on an appropriate choice for e cut , so that the modes at the cutoff have an average 
occupation of order unity. This choice means that all the modes in C are appreciably occupied, 
justifying the use of a classical field for these modes. In contrast the I region contains many 
sparsely occupied modes that are particle-like and would be poorly described using a classical 
field approximation. 

The equation of motion for ij)c is the PGPE 



El ESQ 



tf(r) = Vc(r)+Vi(r), 



(27) 




Wc + MWclVc} 



(28) 



where the projection operator 




(29) 
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formalises our basis set restriction of ip c to the C region. The main approximation used to 
arrive at the PGPE is to neglect dynamical couplings to the incoherent region [26]. The C and 
I regions are then treated as independent systems in thermal and diffusive equilibrium. 

An important feature of Eq. (1281 is that it is ergodic, so that the microstates the classical 
field ipc evolves through in time form a sample of the equilibrium microstates, and time- 
averaging can be used to obtain macroscopic equilibrium properties. 

The procedure for finding equilibrium states including detail on using appropriate input 
states is given elsewhere [12711 . Here, we focus on the method necessary for studying time 
correlations in the system. We define the temporal correlation function for the coherent region 
as 

9c{txM) = (VcMi)V>cM 2 )}, (30) 

i.e. we replace the full Bose field operator with the coherent region c-field. Analogues of G\2 
and Gn follow from this, as previously. 

In previous studies [28] it was shown that the contribution to the spatial phase correlation 
function by the incoherent region was negligible, as any contribution was due only to thermal 
atoms and existed over very short range. Here we neglect the incoherent region when 
calculating phase correlations. Contribution from purely thermal atoms are of course already 
in the coherent region, via the higher energy modes included there. 

We approximate the average in Eq. [30] with an ensemble average. Our procedure for 
calculating phase correlations within the PGPE is as follows: 

• Create an ensemble of 40 input states each with different random fluctuations, and evolve 
these within the PGPE 

• At equilibrium, calculate ^c( r 5 0)^c( r , t) for each of the 40 simulations 

• Perform the average over the ensemble of 40 simulations to calculate (ipc( r i 0)^c(r, *)}• 

5. PGPE Results and Analysis 

In this section we present the results of the application of our PGPE formalism to a cold 
Bose gas. We consider a system of 23 Na atoms in the elongated traps with frequencies 
{u x ,u y ,u) z } = 287r{10, 10, l}Hz and {u x ,u y , u z } = 287r{20, 20, l}Hz . We study two 
temperatures at each trap, corresponding to condensate fractions 45% (T/T c w 0.62) and 10% 
(T/T c ps 0.8). As the condensate is responsible for any long-time coherence in the system, we 
keep the condensate number similar between the two trap symmetries, with N cond ps 45000 
for the cooler system, and N cond ps 30000 for the hotter. 

5.1. Density at Equilibrium 

We begin with a general discussion of the properties of the densities of our four different 
systems. A density slice in the y = plane is shown for each in Fig[T|(see caption). 

The density graphs show one instance of the distribution as it evolves in time. The 
three dimensional PGPE is ergodic, and as such, each simulation evolves through the set of 
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Figure 1. Density slice in the y=0 plane for the four systems studied. Each subfigure is 
normalised separately. From left to right, a) Aspect ratio A = 10, T/T c ss 0.8, b) A = 20, 
T/T c 0.8, c) A = 10, T/T c w 0.62, d) A = 20, T/T c m 0.62 

allowable microstates. In this way, a single shot from a PGPE simulation can be likened to a 
single shot from an experiment - i.e., it gives a microstate rather than an equilibrium average. 
This dynamic quality can easily be seen in the four microstates shown in Fig. [Q where the 
density for each is highly irregular. As well as showing the effect of the changing aspect ratio 
on the basic shape of the cloud, these slices also show the increasing fragmentisation of the 
condensate density at high temperatures or aspect ratios [10]. This is perhaps most striking in 
Fig [Lb, where high temperature and high aspect ratio cause the density to form discrete high 
density clumps, as seen in the separated circular features approximately halfway down the 
figure. 

5.2. Autocorrelation 

The autocorrelation is calculated according to Eq. (|30l) . and the results are presented in Fig [2] 
for our four systems (see caption). 

The autocorrelation for all cases comprises of three distinct parts: an initial sharp decay, 
subsequent larger decay, and long-term quasiperiodic behaviour. 

The initial sharp decay occurs on a timescale t <C tu^ 1 for all systems, and represents 
the initial loss of coherence over short timescales by the thermal modes of the system. This is 
due to the significant proportion of the thermal cloud included in the PGPE description. 



Temporal correlations of elongated Bose gases at finite temperature 



9 



0.85 - 



0.95 



0.8 - 



0.9 





0.75 - 



0.7 - 



0.65 - 



X=20, T/T =0.8 



0.01 



0.02 



0.03 



0.04 



0.05 



t(s) 



Figure 2. Autocorrelation function <?n for four systems of study: variables given in Figure. 



The subsequent decay in autocorrelation is due to the changing phase of the condensate 
in each system with respect to time, and it is this feature of the autocorrelation that differs 
so markedly between each system. It can clearly be seen that the proportion of condensate 
present in a system affects how quickly it will decohere, as can be seen by comparing curves 
of different temperature and same trap (see caption). Perhaps more surprising is the marked 
difference that the aspect ratio of the trap makes to the decay in coherence. For simulations 
with similar condensate fraction and number, but different aspect ratios, the tighter trap leads 
to a much more rapid decoherence, as can be seen by comparing simulations of similar 
temperature and different aspect ratio. For the two T/T c ps 0.62 systems (Figure (2), top 
two curves) the A = 10 system has gn = 0.92 at t = 0.02s, while the A = 20 system has 
gn = 0.82. This is a large difference for the very modest changes in aspect ratio we simulate. 
This more rapid phase decoherence is due to the lack of a spatially extensive condensate in 
the A = 20 case causing greater phase fluctuation along the system length. 

Following this decay, there is quasiperiodic behaviour in the system, shown as a period 
of revival and further decays of the autocorrelation function, on the timescale t > cu^ 1 . This 
is consistent with the prediction from our Bogoliubov anaylsis results. 

There can also be seen on both curves a fine structure noise. This occurs at a frequency 
of u) x , and is due to fluctuations in the tight dimensions of the trap. We note that the PGPE 
description includes the effects and interactions of all appreciably occupied modes, and for 
these systems, there are several modes in the tight direction that are occupied. Since we 
are interested only in coherence along the long axis, we made no further analysis of these 
fluctuations. 

In order to make a quantitative comparison of our results with those of Section [3l we 
plot ln(— ln(<7 12 )) against ln(£), and extract the gradient from the line of best fit. Results are 
given in Tabled! both for very short times f C u; 1 and short times t w cu^ 1 . Comparing 
our results, we see that at very short times, all results correspond approximately to 2, i.e. 
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gn oc e~* 2 as expected. For t ~ u^ 1 , however, it is found that the exponent is different for 
each case. From Bogoliubov analysis, the expected result for the exponent of the decay is 
0.5 for the 3D case, and 1 for the ID case. The results for aspect ratio A = 10 trap confirm 
that this system is still well within the 3D regime. For aspect ratio A = 20 however, it can 
be seen that the results are intermediary to the two predictions. This can be explained in 
terms of the relative temperature to tight aspect ratio in our simulations. For aspect ratio 10, 
the temperature required for excitations in the tight trap direction is T fa ^ = 6.7nK, and 
the temperature of the cooler system is ten times this, and shows clear excitation of the axial 
modes. With an aspect ratio 20, T fa ^ = 13nK, and the simulation presented here has 
a temperature 8 times this value. In this case, it can be seen that more of the axial modes 
are 'frozen out', and the fitting of the exponent gives a result intermediary to the ID and 3D 
limits, as expected. 



A 


T/T c 


t < u- 1 


t fa u- 1 


10 


0.62 


2.0 


0.49 


20 


0.62 


1.9 


0.71 


10 


0.8 


2.0 


0.44 


20 


0.8 


1.9 


0.62 



Table 1. Results of fit to autocorrelation function gn, for four systems, and two time limits. 

It can also be seen in the fitted results that the exponent at times t ~ co^ 1 decreases with 
increasing temperature, leading to a faster decay in gn. This is expected as the condensate 
fraction is responsible for long time coherence, and this is diminished at higher temperatures. 

As a further means of analysis, we calculate the frequency spectrum of the system from 
the autocorrelation function as 

/ drdte^(r c (^t)M^0)). (31) 

J — oo 

The corresponding power spectra for the four systems are shown in Figure[3] For A = 10, at 
low temperature we see a set of well defined peaks close to the ID Bogoliubov frequencies 
given in Eq. (l20l) . At higher temperature a number of features is observed: Broadening of 
peaks, shifting of frequencies, possibly splitting of modes and what appears to be a few new 
peaks that do not correspond to a ID Bogoliubov frequency. For A = 20, one mode appears 
to dominate, corresponding to the mode with quantum number n = 2. This corresponds to 
the longitudinal breathing mode of the system, as would be expected. 

6. Conclusions 

The finite temperature Bose gas is a nontrivial system displaying complex properties of phase 
coherence. The transition between coherence in the ID and 3D limits has been little studied, 
and here we seek to determine the system behaviour in this regime. We develop formalism 
based on the PGPE to numerically study these symmetries, and compare this to our results 
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Figure 3. Results of fourier analysis. Top: aspect ratio A = 10, T/T c w 0.8 (black), 
T/T c w 0.55 (grey), prediction from ID Bogoliubov analysis (dotted). Bottom: A = 20, 
T/T c w 0.8 (black), T/T c w 0.55 (grey),prediction from ID Bogoliubov analysis (dotted). 



from Bogoliubov analysis of the ID and spherical limits. It is found that systems with aspect 
ratio 10 and 20 both give exponents for the decay of the correlation time that are intermediate 
to the ID and 3D limits. This can be explained in terms of the finite temperature of the system 
causing excitations to be present in the tight directions of the trap, placing our systems of 
study deeply within the ID to 3D crossover region. 

For experimental studies of the phase evolution of the systems studied here, a 
splitting and recombination method HU would be impractical, since with reasonably high 
temperatures, splitting and recombining in real space is too slow to capture the initial decay 
of phase coherence. In this case the systems would be better probed with a Bragg-type 
measurement [H HI in which atoms are outcoupled using a light pulse, and the momentum- 
dependent power spectrum is directly obtained. 

This project was financially supported by the Swedish Research Council, Veten- 
skapsradet, and the Kempe foundation. This research was conducted using the resources 
of the High Performance Computing Center North (HPC2N). 
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